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ABSTRACT 

Aims. Use the standard fireball model to create virtual populations of gamma-ray burst afterglows and study their luminosity func- 
tions. 

Methods. We randomly vary the parameters of the standard fireball model to create virtual populations of afterglows. We use the 
luminosity of each burst at an observer's time of 1 day to create a luminosity function and compare our results with available obser- 
vational data to assess the internal consistency of the standard fireball model. 

Results. We show that the luminosity functions can be described by a function similar to a log normal distribution with an exponential 
cutoff. The function parameters are frequency dependent but not very dependent on the model parameter distributions used to create 
the virtual populations. Comparison with observations shows that while there is good general agreement with the data, it is difficult to 
explain simultaneously the X-ray and optical data. Possible reasons for this are discussed and the most likely one is that the standard 
fireball model is incomplete and that decoupling of the X-ray and optical emission mechanism may be needed. 
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Since its discovery in 1997 dCosta et alj 19971: Ivan Paradiis ~et all 
119971) . the standard fireball model, or extensions thereof, has 
been use d to model gamma ray burst (GRB) after glow emis- 
sion (e.g. |Panaitescul2005t|J6hannesson et al.ll2006bl) . With ded- 
icated follow-up observations, the number of detected afterglows 
has been steadily increasing. The rapid burst localisation of the 
Swift satellite and its on-board X-ray (XRT) telescope have sub- 
stantially increased the afterglow detection rate and the number 
of detected afterglows is now well over 20CQ This has inspired 
statistical investigations of afterglow properties; among which 
are st udies exploring the luminosity function at a given observer 
time dGendre & Bqerl 120051: iKann et all 120061: Liang & Zhatiil 
2006; Na rdini et alJl2006l 120071) . " 

Motivated by this previous work, we have created theoreti- 
cal luminosity functions of GRB afterglows using the standard 
fireball model. We show that these luminosity functions can be 
described to a high degree of accuracy by an analytical function. 
Furthermore, we show that the parameters of that function are 
for the most parts independent of the chosen model parameter 
distributions, giving a robust estimate of the luminosity function 
of GRB afterglows. We also show that our theoretical luminosity 
function is similar to the observationally determined one, both 
in X-rays and optical wavebands but only if they are considered 
separately. It is difficult to explain the observed luminosity func- 
tions simultaneously in both wavebands. Several possible rea- 
sons for this are discussed, the most likely one being that the 
standard fireball model is incomplete. 

The outline of this letter is as follows: The procedure for 
creating the luminosity functions is described in section 2 and 
we present our results and discuss how they depend on the model 



Table 1. Default values for the upper and lower limit of the pa- 
rameter ranges. All parameters except for p are distributed loga- 
rithmically in the interval. 



Parameter 


Lower limit 


Upper limit 


Eo [erg] 


10 « 


10 M 


w [cm" 3 ] 


0.1 


10 


0o [deg] 


2 


15 


P 


2.1 


2.6 


Se 


0.05 


0.3 


SB 


io- 5 


IO" 3 



1 see http://www.mpe.mpg.de/~jcg/grbgen.html for the most 
recent list 



parameters in section 3. Section 4 is devoted to comparison with 
observations and we conclude the letter in section 5. 



2. A Virtual World of Afterglows 

T he afterglows are created nu merically with the model described 
in |jo~hannesso n et al.1 12006b). It is based on the standard fire- 
ball jet model, in which the energy is injected instantaneously 
into a narrow jet. To create the afterglows, model parameters 
are selected at random fr om pre-defined distributions (see also 
iJohan nesson et al. 2006a). Each parameter range is determined 
so that it represents results found by fitting afterglow observa- 
tions with the standard fireball model or extensions thereof (e.g. 
iPanaitescu & K umar 2001; Panaitescu 20051 fjoriannesson et al.l 
2006b). Table[T]shows the default parameter ranges used in this 
work. The parameter values can in many cases vary by more 
than an order of magnitude from one afterglow to another. All 
model parameters, except p, are therefore varied on a logarith- 
mic scale to get a more realistic parameter disttibution within 
the chosen range. To concentrate on the intrinsic properties of 
the afterglows, we fix the redshift at z = 1 . 

To calculate the luminosity functions, we use the afterglow 
luminosity density at an observer's time of 1 day post-burst. In 
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Fig. 1. Example luminosity functions obtained in our numeri- 
cal calculations at three different frequencies: 2-10 keV (circles, 
solid curve), /?-band (triangles, dashed curve, blue in electronic 
edition) and 8.5 GHz (squares, dotted curve, red in electronic 
edition). The curves are best fits using equation (fl}. For clarity 
we plot L(p(L) x 0.1 for 2-10 keV. Note that the general shape of 
the luminosity function is frequency independent although A, cr 
and Lq are not. 



fact, any time can be used, bu t this choice is well suited for com- 
pariso n with observations ( Li an g~& Zhan gl 120061: iNardini et alj 
2006). One day after the burst is also late enough in the afterglow 
evolution so that both the optical and X-ray emission can be as- 
sumed to originate in the f orward shock of the standard fireball 
model (Nouse ket alj2006h . To reduce statistical errors, each vir- 
tual world consists of 50,000 afterglows, resulting in statistical 
errors of the order of 1% in the calculated luminosity functions. 



3. Luminosity Functions 

Figure Q] shows typical luminosity functions at three different 
observer frequencies. Overlaid on the numerical calculations are 
fits to the results using the function 



<KL) = C\ 



exp - 



ln(L/L Q y 



2cr 2 



exp - 



La 



(1) 



We have chosen this function because of the parabolic resem- 
blance of the plot shown in fig.[U The last exponential factor is 
needed since the luminosity function is not a pure parabola in 
log-log space. Here, Lq is a characteristic luminosity, defining 
an upper exponential cutoff in the function. Both cr and A affect 
the width of the function, controlling its shape, while cr has a 
slightly stronger effect. C is a normalization constant and is de- 
termined by the other 3 parameters. Our results show that 0(L) 
is an accurate description of the numerical luminosity functions, 
giving ax 2 between 50 and 200 in each case for 100 degrees of 
freedom. Note that the function becomes a log-normal distribu- 
tion with an exponential cutoff when A - I. Our results show, 
however, that a value of A « 2 is in the present case more likely, 
varying between the values 1.5 and 2.5, depending on frequency 
and input model parameter distributions (see fig.|3]l. 

To investigate how the model parameter distributions affect 
our results, three different distribution shapes were tested, a uni- 
form one, a Gaussian one and a triangular distribution with a 
maximum in the center of the parameter range. In the case of 



the Gaussian, the standard deviation was set at l/4th the parame- 
ter range and the distribution restricted to the defined range. We 
found that the shape of the model parameter distributions did 
not change our general results, although the exact values of the 
parameters in equation (fl} do depend on it. Since the function 
parameter values change only slightly between the chosen input 
distributions, we decided to use the triangular distribution for all 
our results below. 

In addition to varying the shape of the model parameter dis- 
tributions, we also tested how the parameter limits affect the out- 
comes. In each case, we changed the limits of only one parame- 
ter at a time. To cover all the basic effects, four tests were made: 
first we changed the upper limit keeping the lower limit fixed, 
next the lower limit was varied keeping the upper limit fixed, 
then the width was changed with a fixed center and finally the 
center was changed with a fixed width. The default values for 
the upper and lower limits are shown in table Q] All of the pa- 
rameters, except p, are distributed logarithmically and the width 
is then defined as the ratio between the upper and lower limit. 
Figure [2] shows the most interesting results from our tests. The 
parameters found from fitting the numerically calculated lumi- 
nosity function with equation ([D are plotted against the param- 
eter that is being varied in each case. The general result from 
these tests, is that the function parameters A and cr are not very 
dependent on the limits of the model parameter distributions, 
while Lo can vary by several orders of magnitude. This allows 
us to conclude that the luminosity function can be presented by 
equation ([]]), setting both A ~ 2 and cr ~ 2 and using Lq as the 
only variable. 

We see from fig. [2] that the most significant model parame- 
ter is Eq, the initial energy release. Other model parameters that 
have significant effects on the parameters of <p{L) are Qq, the ini- 
tial opening angle of the jet, and e e and e#, the fractions of en- 
ergy contained in the electron population and the magnetic field, 
respectively. The effects of the external density «o and the elec- 
tron energy distribution index, p, on the function parameters are 
insignificant compared to the other model parameters. 

The general trend when varying the parameter limits is that 
the value of Lq is correlated with the upper limit (lower limit in 
case of 6>o) of the model parameter range, while cr is correlated 
with the width of the parameter distributions. This is shown, re- 
spectively, in the top and middle panel of fig. [2] Note that the 
changes in the parameters of <p{L) depend in many cases strongly 
on the observed frequency. This can, for example, be seen in the 
value of Lq while varying the upper limit of the e e . There are 
no effects in the radio waveband whereas the effects in optical 
and X-ray wavebands are similar and quite significant. Although 
the value of cr depends most strongly on the width of the model 
parameter distribution, it can in some cases also depend on the 
center value of the model parameter distribution. The trend be- 
tween the center value and cr is, however, not as strong as the 
trend between cr and the width. 

The function parameter that shows the least dispersion when 
varying the limits of the model parameter distributions is A. Its 
value is generally clustered around 2 and shows a greater de- 
pendency on frequency than on the model parameters, while 
its frequency dependence is in fact similar to that of cr and Lq. 
Although A shows the least dispersion when varying the model 
parameter ranges, it is the parameter that is most difficult to con- 
strain in the fitting procedure, resulting in the highest 1-sigma 
error. This can be seen in fig. [3] which shows the value of cr as a 
function of A. It is also clear from the figure that there is an anti- 
correlation between the values of cr and A. Pearson's correlation 
test confirms this and results in correlation parameters between 
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Fig. 2. Shown are the values of some of the function parameters found from fitting numerically calculated luminosity functions with 
equation (Q~|i. See fig.[T]for legend. In each case we vary the limits of the model parameter distribution and show only the parameters 
that have the strongest effects. The default parameter limits are given in tableQ] Top panel: The value of Lq while varying the upper 
limit (lower limit in the case of Oq) of the model parameter distribution keeping the lower limit (upper limit in the case of 9q) fixed. 
Shown are results for the parameters Eq, 9o, e e and eg- We multiply Lq by 10 3 for 2-10 keV. Bottom panel: The value of cr while 
varying the width of the model parameter distribution, keeping the center fixed. Shown are results for £0, #0, e e and eg- Since the 
parameters are changed logarithmically, the width is the ratio between the upper and lower limit. 
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Fig. 3. The plot shows that there is an anti-correlation between 
the values of A and <x found from fitting the numerically cal- 
culated luminosity functions with equation £[). The Pearson's 
correlation coefficient varies from -0.5 to -0.8 depending on fre- 
quency, being strongest in the optical range and weakest in radio. 
See fig.Q]for legend. 



-0.5 to -0.8 depending on frequency. The correlation is strongest 
in the optical but weakest in the radio. 

4. Comparison with Observations 

To compare our results with observations we too k the recent lu- 
minosity compilation from Nard ini et al.l d2006) and converted 
it to luminosity functions using the same method as for our nu- 



merical results. This data is suitable for our use, since it is nor- 
malized to an observer's time of 1 day and a redshift of 1. The 
data is displayed in fig. |4] along with selected luminosity func- 
tions from our calculations. Note that our results are not formal 
fits to the observational data, but rather an eye-guided selection 
of one of our calculated luminosity functions. The sparse data 
set does not warrant a formal fit. Selection effects are also not 
taken into account for the sam e reason, but we refer interested 
readers to Nar dini et alj d2007l) which discusses the selection ef- 
fects. The model parameter distribution that best resembles the 
observed X-ray luminosity function is the one with the center of 
the Eq distribution at 10 51 erg but its width at the default value 
of 100. Other parameters were at their default values, given in 
table Q] This energy distribution i s com parable to the one found 
from observations bv lFrail et al.l ([2001 ). The resulting luminos- 
ity function is, however, incompatible with the observed opti- 
cal luminosity function that shows a clear over abundance of 
high luminosity afterglows which can not be explained by our 
results. There is also a sign of bimodality in the observed optical 
luminosity funct i on dLiang & Zha ng 2006; Nar dini et al.| [2006; 
iKann et al.ll2006tlJakobsson et al.l l2006). 

In an attempt to understand the optical luminosity function, 
we also created two afterglow populations with very narrow 
model parameter distributions, the width being a factor of 2 to 
3. These distributions are not compatible with parameter ranges 
found from fitting afterglow observations with the standard fire- 
ball model, where the ranges seem to be much larger. The two 
populations differ only in the value of the center of the Eq dis- 
tribution and the number of afterglows in each population. The 
more luminous population has the center value of £0 at 10 51 ergs 
while the fainter has the center at 2 ■ 10 50 ergs. The number of 
afterglows in the more luminous population is four times higher 
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Fig. 4. Comparison between our numerical calculations (curves) 
and luminosity fun ctions obtained from observational data 
dNardini et al.l l2006) (filled symbols). Two different models are 
shown, one wide (thick curves), assuming the value of the center 
of the Eq distribution is at 10 51 erg and its width at the default 
value of 100. Other parameters are as in tableQ] The other model 
is a composite luminosity function of two different populations, 
each with a much narrower parameter distribution (thin curves). 
Each population has, in this case, a different value for the center 
of the Eo distribution and a different number of afterglows but 
are otherwise the same. See text for further discussion. Two dif- 
ferent frequencies: 2-10 keV (circles, solid curves) and /?-band 
(triangles, dashed curves, blue in electronic edition) are shown. 



than in the fainter one. Note that this ratio may change if obser- 
vational biases are taken into account, as discussed below. These 
populations give two very narrow luminosity functions, which 
when added together agree with the observed optical luminosity 
function. It does not, however, seem to agree with the X-ray data, 
although a slight hint of a high luminosity peak may be seen. Its 
statistical significance is unknown at present. 

There are a few possible explanations for the difference be- 
tween the optical and X-ray luminosity functions. First, there 
may be two populations of afterglows with different parameter 
distributions, as shown by the thin lines in fig. [4] Since this gives 
a good representation of the bimodal optical luminosity function 
and there is a hint of bimodality in the observed X-ray lumi- 
nosity function, as discussed above, this explanation may seem 
to agree better with the observational data than a single popula- 
tion would. However, the narrow luminosity functions required 
for this explanation force the model parameter distributions to 
be much narrower than those found from afterglow modeling. 
Tight constraints between model parameters might, however, re- 
sult in a similarly narrow luminosity functions and still comply 
with afterglow modeling. But even if that is the case, the ratio 
between the two peaks in each observed luminosity function is 
different and a fit with two afterglow populations would there- 
fore be marginally good at best. 

The second possibility is that the standard fireball model 
needs to be refined. Since there is a distinct difference between 
the optical and X-ray luminosity functions, a decoupling of the 
X-ray and optical radiation mechanism may be needed. One ex- 
ample could be that the X-ray and optical radiation originate in 
different parts of the outflow. The bimodality of the optical lu- 
minosity function might then originate in two different popula- 
tions of afterglows. The lack of bimodality in the X-ray lumi- 



nosity function would then be accounted for by a larger disper- 
sion in the radiation power of the X-ray zone. This possibility 
is also supported by recent observati on of GRB afterglows (e.g. 
iPerlev et all 120071; lOates et alj|2007t) that show incompatibility 
with the standard fireball model because of differences between 
the X-ray and optical light curves. 

Finally, the observed optical luminosity function may be in- 
complete. Limitations of current optical follow-up instruments, 
either in sensitivity or response time, lead to an observational 
bias for luminous afterglows. The bimodality is then expected 
to vanish with increasing number of observations. This may al- 
ready have happened w ith the X-ray afterglow observations, as 
iGendre & Boerl d2005h found early on an evidence for bimodal- 
ity in the X-ray luminosity function, which can now barely be 
seen with the more rapid and sensitive observations. If not real, 
the bimodality in the optical luminosity function may require a 
detection of a larger number of afterglows before it disappears, 
because of the redshift dependence of the data. The redshift is 
in most cases acquired from the optical afterglow spectrum and 
therefore requires an optically bright afterglow. This further in- 
creases the bias for optically bright afterg lows without affect- 
ing the X-rays much. The recent study by Nardini et al. (2007) 
has shown that this is an unlikely explanation. The luminosity 
functions are, however, still based on a very small sample. In 
addition, there are also uncertainties involved when doing the 
K-correction and time shift of the observational data. Hence, we 
cannot completely rule this explanation out. 

5. Discussion and Conclusions 

In this letter we have shown, with the help of numerical calcu- 
lations, that the luminosity functions of GRB afterglows can be 
described by a rather simple function given by equation {T). It is 
similar to a log normal distribution with an exponential cutoff. 
We have also shown that the shape of the luminosity function is 
not sensitive to the shape of the model parameter distributions. 
The function parameters are also quite robust, A and cr vary by 
less than a factor of 2, in most cases only by a few percent, for 
variations of orders of magnitude in the model parameters. It is 
mainly the upper cutoff, Lq, that is affected by the input param- 
eters. This result can be used to get an approximate single pa- 
rameter representation of the afterglow luminosity function. One 
must, however, bear in mind that cr and A are somewhat depen- 
dent on our assumptions for the model parameter distributions. 
While we have chosen those to be representative of parameter es- 
timates from current afterglow modeling, future modeling may 
require a re-evaluation of the values of A and cr. However, these 
are not expected to change much unless significant changes are 
seen in the model parameter distributions or strong correlations 
are found between some of the model parameters. 

We have also shown that our results are in good general 
agreement with observations. However, incompatibility appar- 
ently exists between the optical and X-ray data, where a bi- 
modality is seen in the optical luminosity function that is not 
apparent in the X-rays. The most likely explanation for this is 
that the standard fireball model is an incomplete description of 
the afterglow emission. A more detailed model or possibly an en- 
tirely different one is required, where a decoupling of the sources 
of X-ray and optical radiation is a natural ingredient. This is sup- 
ported by recent observations, which have shown the afterglow 
behav iour to be more complex than expected (e.g. lNousek et alj 
120061) . 
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